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Rough-wall channel analysis using suboptimal 

control theory 

By O. Floresf, J. Jimenez f AND J. Templeton 


1. Introduction 

The original aim of this work was to shed some light on the physics of turbulence 
over rough walls using large-eddy simulations and the suboptimal-control wall boundary 
conditions introduced by Nicoud et al. (2001). It was hoped that, if that algorithm 
was used to fit the mean velocity profile of the simulations to that of a rough-walled 
channel, instead of to a smooth one, the wall stresses introduced by the control algorithm 
would give some indication of what aspects of rough walls are most responsible for the 
modification of the flow in real turbulence. It was similarly expected that the structure 
of the resulting velocity fluctuations would share some of the characteristics of rough- 
walled flows, thus again suggesting what is intrinsic and what is accidental in the effect 
of geometric wall roughness. 

A secondary goal was to study the effect of ‘unphysical’ boundary conditions on the 
outside flow by observing how a relatively major change of the target velocity profile, 
and therefore presumably of the applied wall stresses, modifies properties such as the 
dominant length scales of the velocity fluctuations away from the wall. 

As will be seen below, this secondary goal grew more important during the course of 
the study, which was carried out during a short summer visit of the first two authors 
to the CTR. It became clear that there are open questions about the way in which 
the control algorithm models the boundary conditions, even for smooth walls, and that 
these questions make the physical interpretation of the results difficult. Considerable 
more work in that area seems to be needed before even relatively advanced large-eddy 
simulations, such as these, can be used to draw conclusions about the physics of wall- 
bounded turbulent flows. 

The numerical method is the same as in Nicoud et al. (2001). The modifications in- 
troduced in the original code are briefly described in §2, but the original paper should 
be consulted for a full description of the algorithm. The results are presented in §3 and 
summarized in §4. The elementary properties of turbulence over rough walls which are 
used in the text have been taken from recent reviews such as Raupach et al. (1991) or 
Jimenez (to appear in 2004). 


2. Simulations 

Turbulent channel flow has been simulated using the code presented by Nicoud et al. 
(2001) with Re T = u T h/v = 1000, where u T is the friction velocity and h is the channel 
half-width. The coordinates used are x, y and z, indicating streamwise, wall-normal and 
spanwise directions. The streamwise and spanwise periodicity lengths are L x = 27t and 
L z = 27t/ 3 for all the cases presented here. Except when explicitly noted to the contrary, 
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Re T 

Ax + 

A y + 

A«+ 

L x /h 

L z /h 

Nx X Ny X Nz 

a+ 

A U+ 

SI 

1000 

196 

61 

65 

2n 

2tv/3 

32 x 33 x 32 

6 x 10~ 2 

0 

S2 

1000 

196 

61 

65 

2n 

2tt/3 

32 x 33 x 32 

1.5 x 10~ 2 

0 

R1 

1000 

196 

61 

65 

2n 

2tt/3 

32 x 33 x 32 

6 x 10” 2 

9 


Table 1 . Test cases. The subindices x, y and z indicate the streamwise, wall-normal and span- 
wise directions. The TV’s are the number of grid points in each direction, the A’s are the cor- 
responding grid spaces and the L’s are box lengths. A U + is the roughness function, defined in 
(3.7), and a + is the parameter in equation (2.1). 


all the variables are normalized with u T and h. The superindex + is used for variables 
expressed in wall units. 

The characteristics of the three main simulations used in this paper are summarized 
in table 1. The letter S indicates smooth channels, and R is used for rough ones. The 
reference velocity profile used for SI and S2 is obtained from a direct simulation by del 
Alamo et al. (2003) with Re T = 950. The rough- wall profile is described in the next 
section. 

The LES code uses a standard dynamic Smagorinsky subgrid-scale stress model with- 
out explicit grid filtering. The spatial discretization is a second-order finite-difference 
scheme on a staggered grid, and the time integration is third-order Runge-Kutta with 
an implicit scheme for the wall-normal viscous terms. The flow is driven by a constant 
mean pressure gradient that balances the wall stresses, d x p + = — 1. The overbar stands 
for averaging over wall-parallel planes, while ( ) will be reserved for temporal averaging. 
With the single exception described at the end of §3.1, the time step At is constant. 

The boundary conditions in the streamwise and spanwise directions are periodic, and 
the impermeability condition is imposed at the walls. The two other boundary conditions 
are the total shear stresses r xy and r zy at each point of each wall. They are adjusted at 
each time step to minimize a cost function J , which measures the difference between the 
plane-averaged velocity and a given mean velocity profile U re f. 

J{K, <t>w) = J [(u - U re f) 2 + w' 2 ] dy + a , (2.1) 

where the control variables are vectors of size 2 x N x x N z 

ft U — ( T xy \y=^hi~T xy | y— o) i <j>w = {j zy \y=2hi ~ ly=o)> (2-2) 

and a is a parameter with dimensions [a] = L/U 2 that will be discussed in the next 
section. Primed variables refer to fluctuations with respect to the plane average, so that 
(f>' = <j> — <j>. Note that the global conservation of momentum ensures that (f> satisfies 

® \y=2h -<^} |„=0= A (2-3) 

and 

® \y=2h-{Wy) |„=0=0, (2.4) 

but that the instantaneous mean stresses at each wall may oscillate in time. 

The minimization of J implies an iterative scheme in which an adjoint problem is 
solved in order to compute an estimation for the gradient of the cost function. The 
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Figure 1 . Convergence history of the mean velocity at the first plane off the wall, for a 
smooth- walled channel. , Using the correction (2.5); , without the correction. 

algorithm is suboptimal in the sense that, instead of minimizing the averaged value of J 
over a long period of time, the optimization is computed over each individual time step. 

Several modifications were introduced in the code during preliminary attempts to 
match smooth-wall velocity profiles at different Reynolds numbers. The most important 
was probably the use of the more efficient Brent’s method for the optimization scheme 
(Press et al. 1993), instead of the simpler relaxation used by Nicoud et al. (2001). This 
allowed the stable utilization of the algorithm over a wider range of parameters. 

Another modification involves the recalculation of the mean wall stress. In the original 
code, the value of r™ y + = uf 2 was adjusted at each time step so that the mean velocity 
in the first plane off the wall satisfied the desired logarithmic law, 

u~^~ i 

!°g — A = 0. (2.5) 

The idea was that the mean stress is directly determined by the averaged momentum 
balance, and does not have to be controlled. In each time step the optimum wall stresses 
were computed using the control algorithm, and their mean value was later corrected to 
the u 2 obtained from (2.5). Notice that '«+ = 1 satisfies (2.5) whenever ui =1 is given by 
the desired logarithmic law. While the general argument is sound in a timed-averaged 
sense, this procedure ‘second-guesses’ the control algorithm, since the integrated mo- 
mentum equation is already incorporated in the code, and its effect was found to be 
pernicious. The convergence of the original code was oscillatory and the time history of 
most variables suggested a poorly-damped instability of the control algorithm. Removing 
the correction step suppressed the oscillations and the instability. Two typical conver- 
gence histories, with and without the stress correction, are shown in figure 1. All the 
results below are obtained without using (2.5). 


3. Results 

3.1. The cost function 

A point that needs some discussion is the form of the cost function (2.1). The second 
term on the right-hand side of (2.1) represents the energy of the control variables </>„ and 
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Figure 2. Profiles of (a) the mean velocities and (6) the streamwise velocity fluctuations for a 
smooth-walled channel at Re T fn 1000, using different cost functions, o , case SI; • , case S2; 
, direct simulation (del Alamo et aL2003). 

<j> w , and is introduced to ensure numerical stability. The proportionality coefficient a is 
essentially arbitrary from a physical point of view. 

The value of a used here for the cases SI and R1 is the one given by Nicoucl et al. 
(2001) as being the smallest one that kept their simulation stable. The effect of a larger a 
should be to lower the variance of t w , weakening the strength of the control, and therefore 
presumably degrading the quality of the resulting mean profile. It is nevertheless possible 
to ‘overcontrol’ a system, in the sense that some properties which are not taken into 
account in the cost function may be contaminated by the control itself. It was noted by 
Nicoud et al. (2001) that the velocity fluctuations in the simulated channels are stronger 
than those in direct simulations or in experiments, especially near the wall. Since those 
fluctuations are driven by the wall stresses, which are in turn influenced by a, it is 
conceivable that the reason for the high fluctuations is that a was chosen too small. 

A measure of the strength of the control is the temporal variability of u, compared 
with its ‘expected’ variation. If the velocities are considered to be uncorrelated random 
variables, which is probably not too far from the truth for a coarse LES, the expected 
variance of u would be 

4 = (3.1) 

which has to be compared with the actual measured value 

a 2 = ((u-(u)f). (3.2) 

The ratio a / ap should be of order unity in a natural flow, and its deviation from unity 
is a measure of the strength of the control. 

It is not clear which is the right normalization for a, but the observation in Nicoud et al. 
(2001) that the effect of the control is mainly felt on the velocities near the wall suggests 
that h is not a relevant length scale, and that a more natural normalization would be 
wall units. The value recommended by Nicoud et al. (2001) is then a + k 6 x 104 and 
results in a/ap w2x 10~ 3 . Both the small magnitude of the dimensionless a and of the 
ratio of the standard deviations strongly suggest that the system is overcontrolled. 

A detailed analysis of the effect of a is beyond the scope of this work, and it was 
not attempted, but a simple approximation of the cost function can be used to estimate 
the relevant trends. Assume that the velocity u near the wall is essentially a (random) 
function of the local wall stresses, u = f(r w ), where all the equations in the following 
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argument have to be understood as being applied at a given x and z. We can then 
approximate the cost function as, 

J « 26{u - U ref ) 2 + 2a{r w ') 2 , (3.3) 


where 8 is a measure of the height over which the velocity is modified by the control, 
and u is taken at some representative distance from the wall of order 8. The factor of 
two reflects the contribution from both walls. The minimization of J requires that its 
derivative with respect to all the t w should vanish, giving for each grid point 


j^ = 4 ^-tw _gl + 4a _^ 

dr w N X N Z dr w N X N Z 


= 0 , 


(3.4) 


where we have assumed that t w does not depend of the control. Assuming that f(r w ) 
can be linearized around t w and expressing u in terms of /(r*), we obtain after some 
algebra that 


T 


wf 


[Uref ~ f{T W )] 


8f T /a 
1 + 8f 2 /a 


(3.5) 


where f T is d // d t w evaluated at t w . This is an interesting expression because it captures 
the tendency of t' w to increase with decreasing a, but saturates to some finite value when 
a — > 0. In that limit the cost function becomes u — U re f oc a. 

This has the simple interpretation that, if there is a solution to the unconstrained 
control problem, i.e. if there is a t w such that the first term in (2.1) vanishes, the control 
algorithm approximately identifies it for some small value of a, and any further decrease 
in a does not change that solution appreciably. If such solution does not exist, r ,w and 
v! would diverge as a when a — > 0, and the cost function would saturate to some 
non-zero minimum value. This would appear in our simple model as f 2 = 0. 

We tested this behavior by running case S2, in which a was divided by four with 
respect to the two other cases. The new value of a decreases cr/cxp, as expected, while 
the mean velocity profile improves everywhere (figure 2a), and the intensity of the velocity 
fluctuations near the wall increases (figure 2b). All this is in qualitative agreement with 
the previous analysis. The ratio between the standard deviations of the wall stresses in 
cases S2 and SI is 2:1, which is smaller than the ratio of 4:1 which would correspond to 
an or 1 behavior, and suggests some degree of saturation towards an exact solution. 

To test whether such a solution exists, a second experiment was run with a = 0. The 
velocity fluctuations grew still further, apparently without bound, leading to numerical 
instability for any given constant time step. When the code was modified to run at a 
constant CFL, we were able to trace the growth of the fluctuating velocities in the first 
plane off the wall to values of the order of u ,+ « 600, implying that the suboptimal 
boundary conditions, at least in their present form, are not able to drive the LES to full 
agreement with the desired profile. The errors in figure 2(a) suggest that the difficulty is 
not with the behavior of the near-wall region, but with the center of the channel. It is an 
interesting, although unresolved, question whether this is a limitation of the LES model 
itself, which cannot reproduce the phenomena responsible for the wake component, or of 
the boundary conditions. 


3.2. Rough walls 

The most important effect of roughness on the mean velocity above the buffer zone is 
a constant velocity decrement, A U + (Raupach et al. 1991; Jimenez, to appear 2004). 
A second effect is a shift Ay in the wall-normal coordinate, due to the uncertainty in 
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Figure 3. (a) Mean velocity profiles, (b) Streamwise velocity fluctuations, (c) Wall-normal 

velocity fluctuations. ( d ) Spanwise velocity fluctuations. , smooth wall DNS, Re T = 950; 

, reference velocity profile for the rough case; o , case SI; a ; case Rl. 


the position of the wall. In our case, in which the velocity is only computed at fairly 
large values of y + , the wall-normal shift is negligible and the mean velocity profile can 
be written as 

u + {y + ) = - log(y + ) + 8.5 - -log(fc+) + -H(y/h), (3.6) 

K K K 

or 

U+(y+) = - log(y+) + 5.1 - A U+ + -U(y/h), (3.7) 

K K 

were II is the wake function, and either A U + or kf characterize the effect of roughness 
on the mean velocity profile. 

Using equation (3.7) we construct our rough mean velocity profile for case Rl just 
subtracting a constant A U + = 9 + from the smooth mean velocity profiles in del Alamo 
et al. (2003). This corresponds to a fully-rough flow with an equivalent sand roughness 
k+ « 140. 


3.3. Statistics 

Figure 3 shows the mean velocity and the velocity fluctuations profiles for cases SI and 
Rl. The mean velocities share some characteristics with those obtained by Nicoud et 
al. (2001), even though the latter used as reference profile a logarithmic law without a 
wake function. The error between the mean velocity and U re f is given in figure 4 for 
the three cases computed. It is quite small near the wall, but grows to 0(1) at the 
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Figure 4. Error in the mean velocity profiles with respect to their references, o , case SI; • , 

case S2; a 5 case Rl. 


center of the channel. Between the second and the third grid point there is a weak 
discontinuity that changes the intercept of the logarithmic law. Its origin is probably the 
spatial discretization, because it appears at the same grid point in Nicoucl et al. (2001) 
even if their resolution is coarser than ours by a factor of four. 

The root-mean-squared velocity fluctuations are presented in figures 3(&), 3(c) and 
3(d). The streamwise and spanwise fluctuations near the wall are higher in all the large- 
eddy simulations than in the direct simulation, although the agreement is quite good in 
the core region for the streamwise direction. The wall-normal and spanwise fluctuations 
are slightly underpredicted by the suboptimal code in the core region. The high values 
for the streamwise velocity fluctuations near the wall and the velocity fluctuations in 
the core region agree with the results of Nicoud et al. (2001), but the exceptionally high 
values of the spanwise velocity fluctuations near the wall shown in figure 3(d) do not. 

The lower values of the velocity fluctuations near the wall in case Rl with respect to 
SI are consistent with the effect of roughness on physical turbulence. This is often inter- 
preted as the interference of the roughness elements with the smooth-wall self-sustaining 
turbulence cycle (Jimenez & Moin 1991), as explained by Jimenez (to appear in 2004). 
This interpretation is unlikely in the present case because the resolution is too coarse 
in all three directions to capture the scales associated with the regeneration cycle, and 
the most likely explanation is that the decrease of the fluctuations is due to the lower 
turbulence production near the wall due to the weaker velocity gradients across the first 
grid element of the rough reference profile. 

The wall stresses r xy and r™ y provided by the control algorithm are studied using their 
probability density functions (p.d.f.), which are shown in figure 5. They are compared 
at y + = 60 with the p.d.f. ’s of the shear stresses in a DNS channel with Re T = 550 (del 
Alamo et al. 2003), after box-averaging the latter to the LES grid (A x xA z = 196 + x65 + ). 
Although it was not possible to post-process the shear stresses from the Re T = 950 in time 
for this report, there is probably little differences between them an those at Re T = 550, 
because this part of the flow scales approximately in wall units. 

It is clear from figure 5(a) that the T xy in the suboptimal simulations are different from 
the DNS one, especially in that the former are not able to reproduce the asymmetry of 
the p.d.f. and try to reproduce the mean value of the shear stresses by displacing their 
modes to positive values. As expected, the standard deviations behave in the same way 
as the velocity fluctuations, Rl < SI < S2. Figure 5(c) and 5(d) show the p.d.f. ’s of the 
three cases when the stresses are normalized with theirs mean value (r) and standard 
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Figure 5. Probability density functions (p.d.f.) of the shear stresses at the wall. The shear 
stresses in figures (a) and ( b ) are expressed in wall units, while in figures (c) and (d) they are 

normalized as in equation (3.8). , DNS (see text for details); , Gaussian; o , case SI; 

• , case S2; a ? case Rl. 


deviation 


t n = 


t- (t) 


(3.8) 


The collapse of the three curves on the Gaussian distribution, and the differences with 
the DNS results, suggest that there is little physical information on the shear stresses 
given by the control algorithm. 

Another way of characterizing the shear stresses given by the control is the study of 
their spectral distribution. We can define the spectrum E tJ 


Eij ( k x ,k z ) = (t ijTij), (3.9) 

where Tij{k x , k z ) are the Fourier coefficients of the two-dimensional Fourier transform 
of the wall stresses t-j. The asterisk * indicates complex conjugation and k x ,k z are the 
wave numbers. 

Figure 6 shows these premultiplied spectra for the cases SI and Rl, as functions of 
the wavelengths X x = 2tt / k x and \ z = 2n/k z . Because each spectrum is normalized with 
its maximum, it only contains information about the wavelengths. There is very little 
differences between the two cases. 

It is important to notice that the most energetic modes are located in the large- 
wavelength end of the spectra, specially as regards their widths. The ‘infinitely wide’ 
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Figure 6. Two-dimensional premultiplied spectra of the stresses at the wall as functions of the 
streamwise and spanwise wavelengths, (a) E xy , ( b ) E zy . Shaded contours are from the SI case 
and lines are from Rl. Each spectrum is normalized with its maximum. 


(a, 0) modes for the r xy spectrum sums up to 26% of the energy in the smooth case and 
to 35% in the rough one, but the energy in the ‘infinitely long’ (0,/3) modes is below 1% 
of the total in both cases. On the other hand, for the t™ spectrum the energy contained 
in the (cr, 0) modes is negligible, while the energy contained in the (0, /?) modes sums 7% 
of the total energy in the smooth case and 9% in the rough case. That suggests that the 
spectrum of r™ y is wider, but not too much longer than the simulation box, while the r™ y 
spectrum almost fits in it. The behavior of this quantity in the DNS is not known. 

It is interesting to analyze how the control variable characteristics affect flow variables 
such as the velocities. In order to do that, we have studied the energy spectrum of the 
velocity components, E uu , E vv and E WW1 defined as 

Euu(^xi ^z) — (U J kx,kz'U'kx,kz)T (3.10) 

The corresponding one-dimensional spectra are obtained adding equation (3.10) over a 
certain direction in Fourier space. 

Figures 7(a) and 7(b) show the one-dimensional premultiplied spectrum of the stream- 
wise velocity components at y + = 30, which is the first grid point in the mesh. Although 
the gradient of this variable is determined directly by the control, the agreement between 
the suboptimal cases and the DNS results is reasonable good, at least as much as can be 
expected from the coarse grid being used. The same happens for the wall-normal velocity 
component (not shown here), but not for the spanwise component, whose spectrum is 
shown in figures 7(c) and 7(d). The spanwise velocity has a very energetic mode which 
spans the full width of the box and almost all of its length. When the instantaneous 
velocity fields are examined, we observe fairly regular diagonal bands of positive or neg- 
ative spanwise velocity. This is reminiscent of the instabilities observed by Jimenez et al. 
(2001) in a channel with a porous wall, where they were traced to a coupling between 
a weakly-damped mode of the impermeable channel and the porous-wall condition. In 
that case the result was a series of spanwise rollers spanning the full height and width 
of the channel, and it is conceivable that a similar coupling might result in the present 
structures. 

The spurious perturbation of the spanwise velocity disappears as we move a few grid 
points away from the wall, as can be observed in the spectra in figures 8. This relatively 
fast relaxation away from the wall of the defects of the wall layer, which recalls the similar 
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Figure 7. One-dimensional premultiplied velocities spectra at y + = 30, normalized in wall 
units and plotted as functions of the wavelengths \ x (a, c) and A z ( b , d). (a) and (6), streamwise 
component; (c) and (d) spanwise component. , DNS; o , case SI; • , case S2; a ; case Rl. 


relaxation of the fluctuation intensities in figures 2(b) and 3(6), suggests that the effect 
of the wall is relatively local to small values of y, and that the main role of the boundary 
conditions is to provide a correct intercept for the mean velocity. 


4. Conclusions 

In this work we have simulated three channel flows using the suboptimal-control code 
developed by Nicoud et al. (2001), after some modifications to improve its performance. 
We have seen that the suboptimal control code is not able to reproduce the wake compo- 
nent of the mean velocity profile of real channels. The magnitude of the error depends on 
the parameter a which weights the energy of the control in the cost function. Decreasing 
this parameter decreases the error in the wake, but degrades the agreement of the velocity 
fluctuation intensities near the wall. A simplified analysis of two numerical experiments 
with different values of a suggests that the problem with the velocity profile near the 
channel centerline is intrinsic to the simulation procedure, and not a consequence of the 
boundary condition algorithm. This result, together with the reorganization of the flow 
structures away from the wall, raises the question of whether it is possible to control the 
whole flow just by acting on the wall. 

There are also open questions about the effects of a on the mean velocity profile and 







3 


Rough-wall suboptimal control simulation 


423 





Figure 8. One-dimensional premultiplied velocities spectra at y + = 150, normalized in wall 
units and plotted as functions of the wavelengths \ x (a, c) and A z ( b , d). (a) and (6), streamwise 
component; (c) and (d) spanwise component. , DNS; o , case SI; • , case S2; a ; case Rl. 


on the root mean squared velocity fluctuations, where the influence of this parameter 
seems to be more important. More work is needed in that direction. 

We have shown that the suboptimal code is able to match a synthetic profile for a 
rough turbulent channel, and that the trend of the changes in the root-mean-squared 
velocity fluctuation profiles agrees with the expected results. It is not clear whether this 
agreement is due to a change in the physics of the wall, or just a secondary effect of a lower 
mean velocity gradient. However, the information gathered from the spectra and from 
the p.d.f.’s point to the latter explanation, as there are no clear differences between the 
structures near the wall in the smooth and the rough cases, except for the intensities. The 
structure of the wall stresses introduced by the control algorithm bear little resemblance 
to the corresponding physical quantities. 

There are in all the present simulations a spurious organization of the near-wall span- 
wise velocity into large diagonal modes, which is most likely due to some instability of 
the control procedure, but which also disappears a few grid points away from the wall. 
Its analysis also requires further work. 
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